/*  sdfasinb.cpp  by K.Tsuru May 22 2015 */
// this is a reference function since ver. 2.21
/***********************************************************************
SDouble class with binary splitting.
It provides arcsin(x) using Newton's method i.e. by solving the equation
f(x) = sin(x) - a = 0 (x = arcsin(a)).
It does iteration
x1 = x0 + delta
where delta = (a - sin(x0))/ cos(x0)
sin(x0) and cos(x0) can be obtained at the same time by CosSinBS(x0).

This is not so fast. "Asin(x)" is faster a few times.
***********************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

static const char* const func = "AsinBSNW";
//static const SDouble ONE(1.0);

SDouble AsinBSNW(const SDouble& a) {
    int d1 = DDCompare(a, ONE);
    if( d1 > 0 ) a.SetError(a.DOMAIN_ERR, func, 3150);  // |x| > 1
    if(a == 0.0) return 0.0;    // x = 0
    if(d1 == 0) return a.Sign(3105) > 0 ? MPi2() : -MPi2(); // |x| = 1 return pi/2 or -pi/2

    // see "sdfrsqrt.cpp"
    RealSize C;
    uint max_sz = a.MaxSize();
    int itrmax = howpow2( (DFIGURES*max_sz)/DOUBLE_FIG+1 ) + 6;
    int count = 0;
    uint ef = (DOUBLE_FIG*2u)/DFIGURES, fig = C.EffFigures();
    bool fullPrec = false;  //calclation is done in full precision or not
    if(ef > fig) ef = fig;
    ef = ceilpow2(ef); // ef = 2^n
    SDouble A( Dabs(a) );   // A = |a|
    double AD = doubleD(A);
    SDouble x, y, delta, sinVal, cosVal;
    x = asin(AD);

    do {
        if((ef = C.SetEffFig(ef)) >= fig) fullPrec = true;
        CosSinBS(x, cosVal, sinVal);
        y=A-sinVal;
        delta = y/cosVal;
        x += delta;
        ef *= 2;
        if(ef >= fig) ef = fig;
        C.SetEffFig(0);
        count++;
    } while( ( !delta.IsMLT(ONE) && (count < itrmax) ) || !fullPrec);
    if(count >= itrmax) x.SetError(x.FATAL, func, 34010);
    x.iterationCount = count;
    return a.Sign() > 0 ? x : -x;
}